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PREFACE 


This  report  presents  a  new  analytic  approach  to  optical  pulse 
propagation  through  a  multiple  scattering  medium.  Such  a  model  is 
needed  to  resolve  propagation  issues  in  the  Blue-Green  Optical  Com¬ 
munication  Program.  The  analytic  model  is  capable  of  (1)  duplicating 
highly  aggregated  Monte  Carlo  computations  at  significantly  lower 
cost,  and  (2)  providing  detailed  computations  (e.g.,  involving  finite- 
dimension  receivers)  that  are  impractical  to  simulate. 

The  effort  was  jointly  funded  by  the  Defense  Advanced  Research 
Projects  Agency,  the  Office  of  Naval  Research,  and  the  Naval  Electronics 
System  Command.  The  technical  advisor  was  the  Naval  Ocean  Systems 
Center.  The  report  should  be  immediately  useful  for  verifying  or  up¬ 
dating  modular  expressions  in  the  current  Navy  model  for  single-pulse 
downlink  propagation. 
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ABSTRACT 


V 


This  report  develops  an  analytic  model  for  the  propagation  of 
an  optical  pulse  through  a  multiple  scattering  medium.  Such  a  model 
is  needed  to  investigate  the  effect  of  clouds  on  optical  communica¬ 
tions  from  a  satellite  to  a  submarine.  A  key  initial  result  is  the 
derivation  of  simple  expressions  for  the  first  two  spatial  and  angular 
moments  of  the  radiance  distribution  for  a  narrow  delta-function  source 
immersed  in  an  infinite  scattering  medium.  The  moments  support  a 
diffusion  approximation  for  the  transport  process  in  an  infinite 
plane-parallel  cloud.  First  the  radiance  is  calculated  at  the  cloud 
exit  and  on  a  plane  an  arbitrary  distance  below  the  cloud,  then  power 
collected  by  a  finite  receiver  located  on  this  plane  is  computed. 

The  model  is  validated  by  comparing  its  results  with  computer  simu¬ 
lation  curve  fits  for  optically  thick  clouds  (i.e.,(T )>  15).  The 
model  is  capable  of  duplicating  nearly  all  the  simulation  results  but 
at  significantly  lower  cost.  Furthermore,  detailed  calculations  im¬ 
practical  to  simulate  are  readily  computed. 
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I.  INTRODUCTION  AND  SUMMARY 


This  report  develops  an  analytic  model  for  the  propagation  of  an 
optical  pulse  through  a  multiple  scattering  medium.  Such  a  model  is 
needed  to  investigate  the  effect  of  clouds  on  satellite-to-submarine 
communications  by  means  of  a  blue-green  laser.  Figure  1  shows  the 
portion  of  the  propagation  path  considered  here  as  well  as  the  key 
physical  and  system  parameters  of  the  environment  and  the  receiver. 

We  consider  both  a  narrow  collimated  and  a  broad  Gaussian  beam 
as  a  source  function  at  the  cloud  entrance.  In  both  cases,  however, 
the  source  is  assumed  to  be  a  delta  function  in  time.  The  narrow  beam 
response  provides  both  a  Green's  function  for  the  broad  beam  input  and 
a  tool  for  comparing  our  results  with  Monte  Carlo  simulations.  The 
cloud  is  assumed  to  be  a  plane-parallel  slab  of  infinite  horizontal 
extent.  Cloud  parameters  are  physical  thickness  Zq,  optical  thick¬ 
ness  T,  scattering  coefficient  k  ,  single-scatter  albedo  A,  average 

s  2 

cosine  of  the  polar  scattering  angle  g  =  (cos  0),  and  (cos  0).  We 
provide  analytic  expressions  for  the  radiance  at  the  cloud  exit  (i.e., 
at  z  =  T  =  kgZ0)  and  the  power  received  by  an  infinite  plane  receiver. 

The  latter  expression  corresponds  to  the  problem  addressed  by  computer 
simulation  experiments. 

Figure  1  shows  that  the  radiance  emerging  from  the  cloud  is  propa¬ 
gated  through  the  atmosphere  to  a  plane  just  above  the  ocean  surface 
(i.e.,  z  =  T  +  d).  The  atmospheric  layer  below  the  cloud  is  character¬ 
ized  by  its  physical  thickness  Z^  and  absorption  coefficient  kg;  the  nor¬ 
malized  distance  d  in  Fig.  1  is  given  by  d  =  k  Z. .  We  calculate  the 

S  X 

power  into  an  on-axis  receiver  (at  z  =  T  +  d)  with  aperture  a  and  polar 
field-of-view  half-angle  ®£ov’  ln  an  actual  communication  link,  the 
receiver  would,  of  course,  be  located  on  a  submarine  some  distance 
below  the  ocean  surface.  However,  we  calculate  the  received  pulse 
above  the  ocean  surface  to  develop  modular  expressions  for  transmission 
and  pulse  spreading  due  to  clouds  and  atmosphere  alone,  and  to  enable 
experimental  verification  of  our  results  with  data  gathered  by  ground- 
based  receivers  [1]  or  scaled  laboratory  experiments. 


The  most  complete  mathematical  description  of  the  multiple  scatter¬ 
ing  problem  is  given  by  the  nonstationary  radiative  transport  equation 
of  Chandrasekhar  [2].  Written  in  dimensionless  variables,  the  transport 
equation  for  the  radiance  I(x,  a;  £)  takes  the  form 

3l(x,  a;  0  .  /• 

- 5f - +  *  V  =  _I  +  in  I  I(-’  a’ ;  ?)p(a* ;  a)  da,Ct* 

“4tt 

+  SQ(x,  a;  O  ,  (1) 

where  =  k£ct — t  signifying  time  and  kg  the  extinction  coefficient; 
x  =  k£X  is  the  normalized  cartesian  spatial  coordinate;  a  =  (0,  4>) 
denotes  an  ordered  pair  of  angular  direction;  and  e^  is  the  unit 
vector,  and  du)^  the  solid-angle  element,  in  direction  a.  In  addition, 
p(aT;  a)  is  the  scalar  phase  function,  Sq  is  the  source  function,  and 
X  is  the  single-scatter  albedo. 

Unfortunately,  the  nonstationary  transport  equation  has  not  been 
solved  in  a  form  suitable  for  computation.  Hence  two  procedures  have 
been  widely  used  to  obtain  quantitative  results — solving  Tr_,  (1)  under 
the  small-angle  assumption  [3],  and  simulating  the  scattering  process 
by  computerized  Monte  Carlo  methods  [4]. 

Although  the  scalar  phase  function  is  indeed  highly  peaked  in  the 
forward  direction  for  aerosols  whose  mean  radius  is  much  larger  than 
the  optical  wavelength,  the  small-angle  assumption  may  be  invalid  for 
clouds  thicker  than  several  mean  free  paths.  In  principle,  Monte 
Carlo  simulations  can  provide  the  required  numerical  results  for 
studies  of  optical  communication  systems.  In  practice,  however,  the 
cost  of  simulations  usually  limits  calculations  to  highly  aggregated 
quantities  (e.g.,  optical  pulses  collected  over  infinite  receiving 
planes  [4]).  More  important,  Monte  Carlo  curve  fits  usually  do  not 
illuminate  the  underlying  physical  phenomena  of  interest. 

To  overcome  such  difficulties,  we  begin  by  deriving  the  first 
two  spatial  and  angular  moments  of  the  radiance  function  for  a  narrow, 
collimated,  delta-function  source  immersed  in  an  infinite  medium  (Sec. 
II) .  The  moments  are  expressed  as  simple  closed-form  algebraic  formulae 
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derived  with  no  approximations  (e.g.,  small-angle  assumption,  restricted 
scalar  phase  function)  and  span  the  regimes  from  forward  scatter  (opti¬ 
cally  thin  media)  to  diffusion  (optically  thick  media) .  Their  asymp¬ 
totic  limits  are  shown  to  correspond  to  the  usual  thin  and  thick  clcud 
approximations  (Sec.  III). 

In  addition  to  clarifying  the  scattering  process,  the  moments  also 
provide  the  nonstationary  drift  and  diffusion  coefficients  for  the  dif¬ 
fusion  equation  approximation  to  the  transport  equation  (1) .  The  full 
radiance  distribution  at  the  cloud  exit  (for  x  ^  15)  is  derived  from 
the  energy  density  solution  to  the  diffusion  equation  for  finite  cloud 
boundary  conditions  coupled  with  the  published  angular  distribution 
for  a  diffusing  particle  emerging  from  a  scattering  medium  (Sec.  IV). 

The  Green's  function  atmospheric  propagator  then  allows  power  pulse 
calculations  for  a  finite  receiver  placed  an  arbitrary  distance  below 
the  cloud  (Sec.  V).  Finally,  the  effect  on  the  blue-green  laser  com¬ 
munication  program  is  discussed  (Sec.  VI). 

We  validate  the  diffusion  approximation  to  the  transport  process 
by  comparing  our  calculations  with  published  simulation  curve  fits 
for  the  functional  form  of  received  power  versus  time,  total  trans¬ 
mission  versus  the  optical  thickness  of  the  cloud,  spatial  and  angular 
spreading  versus  optical  thickness,  and  multipath  time  delay  versus 
optical  thickness.  In  each  case,  the  two  methods  show  excellent  agree¬ 
ment:  the  calculated  curve  generally  overlaps  the  simulation  fit  for 
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II.  MOMENTS  OF  RADIANCE  DISTRIBUTION 


Consider  a  photon  immersed  in  an  infinite,  nonabsorbing,  scatter¬ 
ing  medium  with  collision  geometry  shown  in  Fig.  2.  Path  lengths  between 
scatterings  &  are  independent,  identically  distributed  (i.i.d.)  random 
variables.  The  Si  are  distributed  exponentially  witl  j  parameter  of 

unity  corresponding  to  randomly  (i.e.,  Poisson)  distributed  aerosols 

T 

and  normalization  of  x  =  (x,  y,  z)  to  the  mean  scattering  distance 
kg\  Polar  and  rotational  deflection  angles — 9^  and  <}>n,  respectively — 
are  i.i.d.  and  distributed  according  to  the  scalar  phase  function, 
which  is  assumed  to  be  such  that  0n  and  cf>  are  independent  and  (by 
symmetry)  <{>  is  uniformly  distributed  on  (0,  2ir) . 

The  density  function  f(x,  a;  Q  is  defined  so  that  f  dV  do)  gives 

*  X  ot 

the  probability  that  a  photon  initially  heading  in  the  z  direction  at 

x  =  0  will,  at  time  £  =  k  ct,  be  found  in  the  volume  element  dV  =  / 

—  s  x 

dx  dy  dz  while  heading  within  the  solid-angle  element  dui^  oriented  in 
direction  a  =  (9,  <j>) .  The  radiance  I(x,  a;  0  in  dimensionless  vari¬ 
ables  resulting  from  an  initial,  arbitrarily  thin  pulse  in  space  and 
time  is  then  given  by  I(x,  a;  0  =  Egf(x,  a;  O ,  where  EQ  is  the  energy 
in  the  initial  pulse. 

Without  explicit  knowledge  of  f(x,  a;  0 ,  this  section  derives 
expressions  for  the  following  moments  of  the  scattering  process,  all 
conditioned  on  a  photon  initially  at  x  =  0  and  directed  along  the 
z-axis  in  an  infinite  nonabsorbing  medium. 


•  The  mean  penetration  into  the  medium  at  time  0 


z(0  =  If  zf  dVx  du)a 

v  n 

x 
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•  The  variance  along  the  axis  of  propagation: 


a 

X 


a2z(0  =  /  /  [z  -  i(5)]2f  dVx  da)a 


•  The  transverse  variance: 


a2x(0  =  a*(5)  = 


a 

X 


x  f  dV  da) 
x  a 


•  The  average  cosine  of  the  polar  heading  at  time  £: 


(cos  0(5))  = 


a 

X 


cos  0  f  dV  do) 
x  a 


•  The  mean  of  cos  0(5): 


(cos  0(5))  = 


a 


cos  0  f  dV  da) 


a 


where  the  solid-angle  integrals  are  over  =  4ir  steradians 
and  x(5)  -  y(5)  =  0  by  symmetry. 


The  following  derivations  are  included  for  mathematical  completeness, 
and  may  be  omitted  by  those  most  interested  in  the  physical  significance 
of  the  moments  described  in  Sec.  III. 

The  derivations  proceed  in  two  parts.  First,  the  (conditional) 
moments  of  fn(x,  a;  5) — i.e.,  the  volume-angle  density  restricted  to 
photons  that  have  scattered  exactly  n  times — are  found;  then  the  un¬ 
conditional  moments  are  computed  directly  as,  for  example. 


"■TTTag.- 
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i(5)  -]£}  ^{n(5)  =  n}  5n(C)  , 


n»0 


where  the  probability  of  n  scatterings  at  time  5a  0  is  given  by 


^{n(C)  =n^=^Te^>  n  s  0  . 


(2) 


CONDITIONAL  MEAN  x  (?) 

_ — n 

At  time  £  the  random  position  vector  of  a  photon  initially  at 
x  *  0  traveling  in  z  direction  and  scattered  n  times  is 


'V® 


yn<€) )  *  V®  =  E  V® 


V 

— m 


Zn<^ 


m*=0 


(3) 


where 


v 

-m 


■fh(!)  • 

k*0  X  ' 


(4) 


l  is  the  random  distance  traveled  by  the  photon  between  the  m  and 


m 


m  +  1  scattering,  and  B.  is  s  matrix  giving  the  change  in  heading  after 

til  ^ 

the  k  scattering  event.  Therefore, 


r 


.(;:n 

\0  0  1/ 


and 


(c.os  cf)^  cos  8^  -sin  <J>k 
\  *  (  sin  4>k  cos  0k  cos 
-sin  8,  0 


cos  <{>k  sin  8^\ 
sin  4»k  sin  0^ 
cos  8, 


. . -  ~  •  —  - 
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The  independence  of  the  9^  and  (J)^  imply  the  independence  of  the  B^, 
although  S,  is  not  independent  of  &m,  for  all  m,  m '  s  n,  when  con¬ 
ditioned  on  a  total  of  n  scatterings. 

Taking  the  expectation  of  Eq.  (3),  we  have 


v»  -£ 

m«0 


a  v 
m  -m 


(3*) 


To  evaluate  l  ,  we  must  digress  into  order  statistics.  Given  that 
m 

the  photon  has  been  scattered  n  times  at  time  £,  let  T^,  %2>  •••>  Tn 
be  the  ordered  random  scattering  times.  Since  unordered  scattering 
times  are  i.i.d.  and  uniform  on  [0,  £],  the  ordered  random  scattering 
times  have  the  joint  probability  density  [5] 


fT  T 
t1,t2*  * 

T  ( Ki »  ♦ 

*  * ,Tn  1 

r  \  s  Hi. 

’  n  ^n 

0  S  S  ... 

•  ^n  ^ 

=  0  , 

otherwise  . 

The  scattering  lengths  are  related  to  the  order  statistics  by  the 

formula 

\  ~  Tk+1 

-  T.  ,  k  = 

0 )  1)  « • *  j  n  * 

(5) 

Thus, 

,  /  J 

^n 

/• 

5k+2  ?k+l 

/*  /• 

\ 

"{( d?n* 

/  d5n-:  "• 

’a 

/  *1*1  /  d?k 

0  0 

?k 

A  h  h  \ 

*  J  *h-i"J  4h  J  *h I’Tti 
0  0  0  / 


«  -f  '■**-  V  'r 


Wh  3 


Since  the  are  i.i.d,  it  follows  from  Eq.  (4)  that 


where 


Thus, 


v»  *  B0  B1  B2 


Bo(0|-Bm|0 


/  0  0 

B  -  0  0 

\-<sin  9)  0  (cos 


l  ). 

s  e>/ 


V  = 

m 


0 

0 


(6) 


Letting  v=l-g=l- 


(cos  0)“/ 

•  (cos  0) ,  with  the  help  of  Eq.  (3')  we  finally  obtain 


xn(?)  “  yn(5)  =  0 


(7) 


and 


■I—  1_-  (1  ~  v) 


z  (?)  = 

n  +  1  v 


n+1 


(8) 


CONDITIONAL  SECOND  MOMENT  (x  (O  xT(C)) 

- — - - -  —n  — nvsy/ 

Multiplying  the  random  position  vector  in  Eq.  (3)  by  its  trans 
pose  gives 
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X  (£)  Xx(5) 
— n  — n 


-tev°4Ev*> 


-E  E  v» 


V  V,  . 

— m  — k 


m=0  k=0 


Taking  expectations  and  using  the  independence  of  the  and  yields 


<2n®  >£«»  Um  V^m  \> 


m=0  k=0 


n  m-1 


"E  ^  +  2E  E  U«  \> 


■  m 


m=l  k=0 


From  Eq.  (5)  we  have 


2  2 
r  =  (x  ..  -  T  r 
m  nrfl  m 


l  £,  =  (T  . .  -  T  )  (T.  ,  -  -  T  ) 
m  k  m+1  m  k+1  k 


^k  ^3  ^2  \ 

/  «k-i  -  /"  «2  /  «i)- 


k  ±  i?r 

■  1)  (n  +  2) 


V 


m8ip3?f!%?£$®i fiSP^IP 


l  I  • 


Substituting  the  above  results  into  Eq.  (9)  and  separating  the 
zeroth  term  in  the  first  summation,  we  obtain 


<x,(S)  £(o>  = 


(n  +  1)  (n  4-  2) 


0  0  0 


\  T""'  T 

0  0  0  1  +  >  (v  ;/> 
0  0  1/  ““  ^ 


n  m-1 


EE(«>  • 


m=l  k=0 


From  Eqs.  (A)  and  (6),  we  have 


^m  ^  ^1  *  ”  Bk  Bk+1  ®k+2 


H!  h} 


-<Br-Bk<Bk+iBk+2-*Bm>(;K> 


.  Is 


=  <BX  B^-MojvJ) 


=  •••  Bk<cos  6)m  kMv£> 


(cos  e)®"1*  (v^  v£> 


m~k  /  T\ 

=  8  <\  V  ' 


it 
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The  evaluation  of  (v^  v^)  proceeds  by  taking  expectations  from 


the  center  outward: 


<4$ 


<B1  B2 


/0  0  0\  T  T  m 


0  0  0 


<»1  <»2  -  <\-l  <M  °  «  °)  "k>  Bll>  -  B2>  BI> 


„T\  t-Tv 


(■^(sin^  0)  0  0  \ 

0  |<si„2e>  0  J  bJ_j> 

0  0  (cos  0)/ 


i  -  r 


0  1  -  C 


0  0 
_k 


0  1  +  2C 


where 


3  (cos  0)  -  1 
5  2 


Since 


?  ’ 


Eq.  (12)  becomes 


rtV ,  v  -  »u 


(xn(S)  yn(£))  (xn  (£)  z. 

(y^(?)>  (yn(0  z. 

<*»«>  yn(S»  <4<E 

Thus,  after  combining  terms  in  Eq.  (14),  we  conclude  that 


*A<c; 
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But  by  definition 


B„  - 


Matching  entries  in  the  last  two  expressions  for  B  ,  the  expectation 

H 

of  the  polar  angle  cosine  after  n  scatterings  is  given  by 


(cos  0)n  =  (cos  0)n  =  gn  =  (1  -  v)n  .  (17) 

Similarly , 

T  /°  0  °\  T 

<s„  V  -  <B«  [l  J  °)  V 

(1/2 (sin2  6)n  0 

0  l/2(si-i2  0> 

11 

0  0 

so,  after  matching  terms  with  Eq.  (13), 


ifi' 


vt. 


(cos2  0)  -  4  (1  +  2£n)  -  4  [1  +  2(1  -  W)n]  .  (18) 

no  -3 

UNCONDITIONAL  MOMENTS 

We  calculate  the  unconditional  moments  from  the  conditional 
moments  after  observing  that  the  probability  of  a  photon  undergoing 
exactly  r„  scatterings  at  time  £  is  Poisson-distributed  with  mean  £, 
as  in  Eq.  (2).  Hence  from  Eqs.  (7)  and  (8),  we  obtain 


and 
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n}  in(C 


=  0  , 


n=0 


y(?)  =2  ^{n(0  = 


n)  yn(?) 


=  0  , 


n=0 


;<?)  =5Z  ^{(n(5)  =  n}  *n(?) 


n*0 


e~g  gn  g  1  -  g  -  V)n+1 


n!  n  +  1 


n=0 


..  £5  ; Y'  £n+1  _  \£  [g(l  -  v)]n+1 

v  \Zu  (n  +  1)  !  Zm^  (n  +  1)  ! 

V  n=0 


n=0 


-£ir.«- 
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?(l-v) 


1  -  (es'*  v'  -  1) 


1  -  e 


-v? 


(19) 


(20) 


(21) 


Similarly,  to  calculate  the  unconditional  spatial  variance  a  (5),  we 


combine  Eqs.  (2),  (15),  and  (19)  to  obtain 


r  -A'V^ 


} 
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III.  MOMENT  ILLUSTRATIONS  AND  ASYMPTOTICS 


Here  we  examine  the  physical  significance  of  the  moment  formulae 
as  well  as  their  behavior  for  small  and  large  times  to  allow  compar¬ 
ison  with  published  results  and  to  demonstrate  the  transition  be¬ 
tween  the  diffusion  and  forward  scatter  regimes.  Our  discussion  is 
predicated  on  the  assumption  of  an  arbitrarily  thin,  delta-function 
time  pulse  immersed  in  an  infinite  nonabsorbing  medium.  Section  IV 
extends  our  results  to  a  finite  cloud  with  absorption,  and  Sec.  V 
considers  a  broad  Gaussian  beam  at  the  cloud  entrance  and  an  atmo¬ 
spheric  layer  below  the  cloud. 

DIFFUSION  REGIME 

It  is  convenient  to  consider  a  simple  diffusion  (e.g..  Brownian 
motion)  as  the  limiting  process  obtained  when  a  particle  suffers  an 
arbitrarily  large  number  of  isotropic  scatterings  separated  by  arbi¬ 
trarily  small  path  lengths.  The  moments  of  the  particle  distribution 
are  then  consistent  with  a  stationary  centroid  at  x  =  0,  a  spatial 
variance  along  each  coordinate  that  grows  linearly  with  time,  and  an 
isotropic  heading  independent  of  time.  Although  individual  photon 
deflections  within  a  cloud  are  highly  anisotropic,  after  a  photon 
undergoes  many  scatterings,  its  behavior  may  be  expected  to  resemble 
that  of  a  diffusing  particle.  This  intuitive  notion  is  verified  by 
the  moment  formulae  below. 

Figures  3  and  4  illustrate  the  temporal  behavior  of  the  first  two 
spatial  moments  of  the  scattering  process.  Figure  3  plots  z(£)  com¬ 
puted  from  Eq.  (21)  for  several  values  of  (cos  9).  The  asymptotic 
limit  of  z(£)  is 


lira  z(£) 

£-KO 


1 

1  -  (cos  0) 


(27) 


The  centroid  of  the  photon  density  is  thus  stationary  at  about  z  =  7 
for  £  >  20  when  (cos  0)  =  0.85,  a  representative  value  for  cloud 
aerosols  [1].  Figure  4  plots  Ox(£"/  and  Oz(£)  from  Eqs.  (22)  and  (24), 


♦ 


standard  deviations  plotted  f< 


Stp$gr9>>5&'* v^tbfr^w***  _*^ 
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2 

respectively,  for  <cos  6>  =  <cos  6)  =  0.85.  The  asymptotic  limit  of 
the  corresponding  variances  is  given  by 

^-^•iiH.se)  <28> 

for  large  C*  Equation  (28)  evaluated  for  (cos  9)  =  0  (isotropic 
scattering)  is  the  classical  diffusion  result  for  Brownian  motion  [6]. 
Equations  (27)  and  (28)  also  evince  the  1  -  (cos  8)  factor  suggested 
for  scaling  diffusion  results  [4]. 

To  conclude  that  the  photon  distribution  for  large  times  resembles 
that  of  a  diffusing  particle,  we  must  inspect  the  behavior  of  the  angular 
as  well  as  the  spatial  moments.  From  Eqs.  (25)  and  (26)  in  Sec.  II, 

(cos  8(C))  =  e  (29) 

and 

acos  e(?)  =  <CO£j2  0(S»  -  <=os  0(5»2 

=  f  (1  +  2e”W?)  -  e_2v?  ,  (30) 

where  0(C)  is  the  polar  direction  of  the  photon  at  normalized  time  C> 

Figure  5  plots  (cos  0(0)  and  a  a(C)  for  (cos  0)  =  (cos2  8) 

cos  tt 

=  0.85.  The  asymptotic  limit  of  (cos  0(C))  is  zero,  corresponding  to 
a  diffusion- like  isotropic  polar  heading;  and 


lim  a 


2 

cos  8 


y 


which  is  also  a  diffusion  result.  Thus,  the  asymptotic  limits  of  the 
first  two  spatial  and  angular  moments  of  the  photon  distribution  are 
those  of  a  diffusion  process  in  an  infinite  medium.  Furthermore,  the 
diffusion  limits  are  approximated  after  about  C  =  15,  suggesting  a 
diffusion  (or  thick  cloud)  regime  for  optical  thicknesses  greater  than 
T  -  15. 


j-.-ggy  ..--MrJFS 


Pia.  5— First  two  angular  moments  tor  <cos  0> 
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FORWARD  SCATTER  REGIME 

Since  photon  scattering  in  clouds  is  highly  peaked  in  the  forward 
direction,  photons  deflected  only  a  few  times  maintain  their  initial 
spatial,  angular,  and  temporal  cohesion.  In  cur  normalized  units,  the 
average  number  of  scatterings  at  time  £  is  5  itself  for  randomly  dis¬ 
tributed  aerosols.  For  small  £,  therefore,  we  would  expect  the  moment 
formulae  to  evince  the  behavior  attributed  to  the  multiple  forward 

scatter  (or  thin  cloud)  regime. 

2 

Expressing  from  Eq.  (22)  as  a  power  series  of  £,  we  obtain 


- 


2  w 


3(w  -  v) 


(-l)k(vk~2 


k-2.  jS* 
"  w  >  kT 


k=3 


V 

which  provides  the  interesting  small-time  expression 

a2x(0  a  |  £3  .  (31) 

The  usual  definition  of  mean  square  spread  [3]  results  in 

<r2>  =  <x2>  +  <y2> 

(32) 

for  small  times.  When  written  in  the  dimensioned  variables  R  =  kg^r 

2  2 

and  t  *  (kgc)  £,  with  the  small-angle  assumptions  (cos  0)  as  1  -  (0  ) 
and  etas  Z,  Eq.  (32)  becomes 

<R2)  =  XkeZ3(02>  , 

which  is  identical  to  the  small-angle  scattering  result  (cf.  f3]  and 
[7]). 


- -  -.  --  -  - -  ---  -  - 


|p|^pFgg*> 


“r.w*- 


2  3 

However,  the  dependence  of  (r  )  on  £  in  the  approximate  equation 
(32)  is  valid  for  only  about  £  £  5,  as  shown  in  Fig.  6,  where  the  rms 
beam  spread  V (r^)  is  plotted  using  ^X(C)  from  both  the  exact  equation 
(22)  and  the  approximate  (small  time)  equation  (31) .  We  thus  conclude 
that,  for  calculating  mean  square  spread,  the  small-angle  assumption 
is  questionable  for  optical  thicknesses  greater  than  T  =  5. 

The  moment  formulae  also  allow  for  a  simplified  calculation  of 
the  average  multipath  pulse  delay  for  thin  clouds.  Let  A?  =  Ej  -  T  be 
the  random  delay  between  the  time  an  unscattered  photon  passes  through 
the  plane  z  =  T  (in  normalized  units,  this  time  is  T  itself  since  the 
normalized  velocity  is  unity)  and  the  random  time  £  at  which  the  possibly 
scattered  photon  passes  through  the  plane  z  =  T. 

It  is  reasonable  to  estimate  the  average  value  of  £  as  the  time 
taken  by  the  centroid  of  the  photon  density  along  the  z-axis  to  pass 
through  the  z  =  T  plane.  Thus,  inverting  Eq.  (21)  yields 


5“V‘”<i-r7F> 


« -  - T  * 


an  expression  that  should  be  applicable  for  small  T. 

Figure  7  plots  Eq.  (33)  and  the  thick-cloud  simulation  curve  fit  [4] 
for  A£  versus  t.  The  expression  given  by  Eq.  (33)  appears  to  diverge 
from  the  simulation  prediction  (compared  with  a  true  cloud  exit  but  for 
relatively  large  t)  at  about  x  =  5,  where  backscatter  starts  to  become 
important.  The  divergence  arises  because  Eq.  (33)  for  A?  is  derived 
from  an  infinite  cloud  model,  so  photons  that  would  have  left  a  finite 
cloud  tend  to  increase  the  time  for  the  z-axis  centroid  to  reach  the 
plane  z  =  x.  For  small  x,  however,  nearly  all  photons  are  heading  in 
the  forward  (increasing  z)  direction,  so  the  lack  of  a  true  cloud  exit 
at  x  is  inconsequential. 

3/2 

An  analysis  based  on  a  path-length  integral  along  r  ~  z  — the 
small- angle-derived,  mean  square,  transverse  spread — obtains  the 


NORMALIZED  RMS  BEAM  SPREAD 


ANALYTIC  APPROXIMATION  A-  SIMULATION  FIT  FOR 
FOR  SMALL  r  P  LARGE  r  (BUCHER) 
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expression  [8]  (in  our  notation  with  the  single-scatter  albedo  set 
to  unity) 

M  ■  ts  7$  K1 + *  t(e,2>)3/2  ■  L]  •  <34: 

Expanding  Eq.  (34)  in  T,  and  ignoring  third  and  higher  order  terms, 
we  obtain  the  small  T  expression 

A?  =  0.57  t2<02>  .  (35 

2 

Similarly  expanding  Eq.  (33)  and  setting  v  -  1  -  (cos  9)  ss  (0  /2) 
yields 


M  -  |  xV>  , 

about  a  factor  of  two  less  than  Eq.  (35).  This  discrepancy  is  appar¬ 
ently  due  to  the  path  chosen  to  calculate  A?  in  Eq.  (34)  (i.e., 
r  ~  z3/2)  being  too  long  for  photons  that  have  suffered,  on  average, 
only  T  scatterings.  However,  Eq.  (34)  eventually  (for  T  >  30)  yields 
an  underestimate  of  A?  compared  with  simulation  results  (cf.  [9]) 
because  the  smooth  path  is  shorter  than  the  jagged  diffusion-like 
sample  paths. 

CONDITIONAL  MOMENTS 

Recall  that  intermediate  results  in  the  moment  derivations  of 
Sec.  II  are  expressions  for  the  moments  conditioned  on  exactly  n 
scatterings.  The  conditional  moments  provide  physical  insights  into 
the  scattering  process  that  are  not  apparent  froui  inspection  of  the 
unconditional  moments.  Hence,  Fig.  8  plots  zr(£)  from  Eq.  (8)  and 


0  (n,  O  =  and  0  (n,  Q  =  [<z2(S)>  -  zhoiH 


NUMBER  OF  SCATTERINGS 


from  Eqs.  (15)  and  (16),  respe  Cively.  Although  Fig.  8  is  for  E,  =  100, 

observe  that  z  (0  ,  O  (n,  0 ,  and  0  (n,  0  are  all  proportional  to  £ — 
n  x  z 

so,  when  considered  as  functions  of  n  for  fixed  0  their  form  is  in¬ 
variant  with  0 

The  behavior  of  z^(0  in  Fig.  8  is  intuitive:  the  higher  the 
scattering  order  n,  the  closer  the  centroid  of  the  conditional  dis¬ 
tribution  lies  to  the  origin  (indeed,  lim  z  (0  =  0).  That  is,  on 

n -*»  n 

average,  the  more  a  photon  is  scattered,  the  farther  it  lags  behind 
its  unscattered  brethren. 

The  plots  of  0  (n,  0  and  0  (n,  0  display  maxima  at  about  n  =  10 

Z  X 

and  n  **  15,  respectively.  The  values  of  n  at  which  these  maxima  occur 

are  sensitive  to  (cos  9)  but  invariant  with  0  Furthermore, 

lim  0  (n,  0  =  lim  0  (n,  0  =0.  Thus,  at  any  time  0  photons 
n-*oo  z  n-*»  x 

scattered  many  times  tend  to  accumulate  at  the  origin,  their  disper¬ 
sion  decreasing  with  order  n.  Such  clustering  occurs  because,  at  a 
given  time,  all  photons  have  the  same  path  length  but  the  relatively 
jagged  (diffusion-like)  paths  of  those  that  have  been  scattered  many 
times  will  have,  on  average,  kept  them  closer  to  the  origin.  That 
effect  would  be  most  pronounced  for  isotropic  scattering  (i.e., 

(cos  9)  «*  0)  since  it  is  identical  to  diffusion  in  the  limit  of  a 
decreasing  mean  free  path  between  scatterings. 


v,„,/ 


i 


3 

| 

! 


-35- 

IV.  DIFFUSION  APPROXIMATION  FOR  MULTIPLE  SCATTERING 
IN  AN  INFINITE  CLOUD  LAYER 


The  moments  of  the  radiance  distribution  and  the  ultimately  dif¬ 
fusive  behavior  of  scattered  photons  are  now  used  to  develop  an  expres¬ 
sion  for  the  irradiance  at  the  exit  of  an  infinite  plane-parallel  cloud. 
The  irradiance  allows  us  to  calculate  the  power  collected  over  an  in¬ 
finite  plane  receiver  at  the  cloud  exit,  then  compare  our  results  with 
Monte  Carlo  simulations  of  the  identical  problem.  The  calculated  and 
simulated  power  pulses  strongly  agree  for  t  =  30  and  are  virtually 
indistinguishable  for  T  =  80. 


5 


£ 


DIFFUSION  EQUATION 

The  three-dimensional  diffusion  equation  for  a  particle  with 

homogeneous  diffusion  coefficient  D(0  =  (D  (0 ,  D  (0 ,  D  (0)T, 

x  ip  y  z 

drift  coefficient  a(0  =  (a  (5) ,  a  (0 ,  a  (0)  ,  and  absorption  co- 

x  y  z 

efficient  K  is  [10]* 


3 P(x,  0 
3? 


K 


p(x>  0  , 


(36) 


where  p(x,  0  is  the  probability  volume  density  for  the  particle  at 
time  0  The  solution  of  Eq.  (36)  for  a  particle  initially  at  the 
origin  (p(x,  0)  =  6(x))  in  an  infinite  homogeneous  medium  is  given 
by  [10] 


p(x,  0  =  (2t r)“3/2  A"1  exp  j-[ic£  +  (2A)  1  (x  -  y)T(x  -  y)]  }  ,  (37) 


where 


and 
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Ji  -  f  a(e’)  dr 

*'0 


A  =  /*  Dx(D  d?’  /*  Dy(£*)  d?’  /*  Dz(5’)  dr  . 

•'O  •'o  ■'o 


Equation  (37)  is  simply  the  product  of  three  independent  Gaussian 


densities  and  an  absorption  term  e 


-k5 


Hence  we  may  write 


p(x,  O  =  e-K^  px(x,  O  py(y,  r  Pz(z,  O  , 


where  p  ,  p  ,  and  p  are  one-dimensional  Gaussian  densities.  Since 
x  y  z 

a  Gaussian  density  is  determined  by  its  first  two  moments  and  the 
moments  of  the  transport  process  in  an  infinite  medium  were  derived 
in  Sec.  II,  we  conclude  that 


3*(5>  / 

a(Q  ~  =  (°»  °> 


Mil 
3C  , 


(38) 


and 


Rtt) 


'3c*(0 

~~W~: 


H  ’ 


da2z(0 


(39) 


where  z(?),  0^(5)  =  0^(?),  and  are  given  by  Eqs.  (21),  (22),  and 

X  y  2 

(24),  respectively. 

-  2  2 

Recall  that  Sec.  Ill  shows  that  z(Q,  0x(£),  and  0z(£)  asymptoti¬ 
cally  approach  classical  diffusion  limits.  Although  the  nonstationary 
transport  equation  (1)  has  been  approximated  by  the  diffusion  equation 
(36)  in  prior  treatments  (cf.  [11]  and  [12]),  the  drift  and  diffusion 
coefficients  have  been  time-invariant — generally  set  equal  to  their 


<m 


;  Sr  § 


asymptotic  values  or  otherwise  approximated.  Armed  with  the  exact 
closed-form  expressions  derived  in  Sec.  II,  we  will  strive  to  extend 
the  region  of  validity  for  the  diffusion  approximation  by  employing  the 
"correct"  time-varying  coefficients.  In  the  simplest  sense,  Eqs.  (37) 
through  (39)  together  constitute  the  correct  Green's  function  for  the 
photon  density  of  the  infinite-medium  transport  problem  because  they 
provide  the  correct  spatial  moments  for  all  £  >  0.  However,  the  dif¬ 
fusion  approximation  used  to  obtain  that  solution  forces  the  density 
p(x,  £)  into  a  Gaussian  form  that  is  valid  only  after  the  onset  of 
the  diffusion  regime  (£  >  15). 

BOUNDARY  CONDITIONS 

Consider  a  photon  that  traverses  the  boundary  between  a  scattering 
and  a  nonscattering  medium.  Since  it  can  never  return  to  the  scattering 
medium,  its  behavior  inside  the  cloud  may  be  likened  to  that  of  a  dif¬ 
fusing  particle  attaining  an  absorbing  wall.  Diffusion  theory  provides 
the  absorbing-wall  boundary  condition  p(x,  £)  =  0  for  all  x  along  the 
wall  [10].  However,  close  to  boundaries,  the  diffusion  approximation 
to  the  transport  process  is  known  to  be  questionable. 

By  direct  solution  of  Milne's  equation  for  a  simpler,  related 
problem,  Morse  and  Feshbach  [11]  have  shown  that  the  appropriate  con¬ 
dition  for  a  diffusing  particle  used  to  approximate  photon  transport 
across  an  exit  boundary  is  given  by 


Pz(z) 


z=0  or  z=T 


where,  as  shown  in  Fig.  1,  we  have  modeled  our  cloud  as  a  slab  of 
infinite  horizontal  extent  and  (optical)  thickness  t  =  k^^ — Zq  being 
the  physical  thickness  of  the  cloud.  Since  Milne's  problem  considers 
a  stationary  plane  wave  incident  on  the  scattering  medium,  Eq.  (40) 
is  the  steady-state  boundary  condition  for  a  three-dimensional  problem 
that  reduces  to  one  dimension  by  symmetry.  Furthermore,  Eq.  (40)  holds 
for  only  isotropic  scattering,  suggesting  that  the  usual  v  =  1  -  (cos  8) 


ill#* 
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scaling  for  anisotropic  scattering  [11-13]  should  be  applied  to  obtain 
the  boundary  condition 


Pz(z)  = 


0.71  8pz 
v  3z 


z=0  or  z=T 


(41) 


which  we  use  below. 

The  difficulty  in  applying  Eq.  (41)  lies  in  specification  of  both 
the  density  and  its  gradient  at  the  boundary.  Adopting  the  technique 
of  Morse  and  Feshbach  [11,12],  we  extend  the  scattering  medium  a  dis¬ 
tance  0.71/v  beyond  each  boundary  and  set  P2(z,  O  =  0  at  the  boundary 
of  the  extended  medium,  as  represented  in  Fig.  9.  Using  the  extended 
medium  with  the  diffusion  approximation  for  Milne's  steady-state  prob¬ 
lem  results  in  an  error  of  no  more  than  a  0.3  percent  for  pz  in  the 
scattering  medium,  compared  with  the  exact  solution  [11* . 

We  now  solve  Eq.  (36)  with  the  coefficients  from  Ecs.  (38)  and 
(39) ,  initial  condition 

p(x,  K)  =  <5(x)  6(y)  6(s)  ,  5-0  (42) 

and  boundary  conditions 

p(x,  0=0,  z  - —  and  z  =  T  +  -- —  (43) 

and 


p(x,  0=0,  x  =  y  =  ± 00  .  (44) 

The  separability  of  the  initial  and  boundary  conditions  for  Eq.  (36) 
implies  that  the  solution  is  of  the  form 

p(x,  O  =  e"K^  Px(x,  O  Py(y>  O  Pz(z,  O  ,  (45) 

where  p^,  p  ,  and  pz  are  the  solutions  to  three  one-dimensional  dif¬ 
fusion  equations  with  the  corresponding  one-dimensional  boundary 
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Fig.  9— Representation  of  photon  density  pz(z,  £) ,  £  fixed,  using  diffusion  theory 
approximation  with  extended  cloud  boundary  condition 
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conditions  and  without  an  absorption  term  [14].  Tims,  p  and  p  remain 

2  2  x  y 

Gaussian  with  mean  zero  and  variance  a  (£)=  o  (£)  given  by  Eq.  (22), 

x  y 

since  our  finite-thickness  cloud  is  still  of  infinite  horizontal  extent. 
Hence 


and 


Pv(x,  O  - - —  exp 

a 

x 


2 

x 


2°k  «> 


P  (y,  £)  =  - - -  exp 

y  a  V^if 

y 


2 


For  Pz(z»  5) »  we  seek  the  solution  of 


D,(0  9  Pz 

2  9z2 


-  az(?) 


9z 


subject  to 


(46) 


(47) 


(48) 


Pz(z»  S)  =  5(z)  ,  t,  =  0 


and 


_  ,  n  0.71  ,  _  .  0.71 

p  (z,  tj  =0  ,  z  - - and  z  =  T  H - . 

2  V  V 

The  method  of  images  obtains  the  solution  [14] 


(49) 
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which  is  a  rapidly  converging  series  for  all  z  and  £  of  interest. 
Thus,  the  solution  of  Eq.  (36)  with  coefficients  given  by  Eqs.  (38) 
and  (39),  initial  condition  by  Eq.  (42),  and  boundary  conditions  by 
Eqs.  (43)  and  (44)  is 

p(x,  £)  =  e"^  px(x,  5)  Py(y,  O  pz(z,  O  , 

where  p^,  p^,  and  Pz  are  given  by  Eqs.  (46),  (47),  and  (49),  respec¬ 
tively.  The  absorption  coefficient  K,  normalized  to  the  mean  free 
scattering  length  kg  ,  is 


K 


(50) 


where  is  the  unnormalized  absorption  coefficient  and  X  is  the  single¬ 
scatter  albedo. 


IRRADIANCE  AT  CLOUD  EXIT 

Consider  che  probability  that  an  unabsorbed  photon  initially  at 
x  =  0  has  left  the  scattering  medium  by  time  £.  This  probability  is 
given  in  the  diffusion  approximation  by  th  cumulative  probability  that 
a  diffusing  particle  remains  unabsorbed  in  che  medium  until  time 
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We  can  thus  use  Eq.  (48)  to  obtain  the  exit  time  (or  "first  passage") 
density 


85 


?{z(Z)  6  [0,  t]>  =  1 


■L 

1 

£ 


•T  3p„ 


35 


dz 


D  (5)  3  P-  3p 

— - t-  -  a  (?)  a-2, 

2  9z^  z  3z 


dz 


-[ 


»  (C)  3p 

az®  ^a-^2-3rJ  -K«)»2 

Z=T 


1  r  ^  V«  3pzl 
j  -Laz(°pz'T -irj  . 

Jz=T  u  Jz=0 


The  first  term  in  the  last  expression  is  the  first-passage  probability 
density  for  photons  leaving  the  cloud  (bottom)  at  z  =  T.  That  density 
is  recognizable  as  the  usual  diffusion  flux 


D 

V  p 
2  z  z 


augmented  by  a  drift  term  a  p  =  e  p  [from  Eqs.  (21)  and  (38)], 

z  z  z 

which  quickly  expires  for  optically  thick  clouds.  Denoting  the  first- 
passage  density  at  the  cloud  bottom  by  A(t,  £) ,  we  have 


A(t 


-»5  V«  s» 

e  p - ; 


fa 

3Z  JZ=T 
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a  (?)  v^2tt  , 
z'  k=-°° 


exp 
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20^0 


r  D  (5) 

^  +  -V-  h 

20‘(5) 


-  exp 
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2  ^ 


2°2z(5). 


c  D  (O 

+  -V-  e2 

2<(5)  . 


(51) 
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where 


\  =  -2(C)  +  T(1  +  2k)  +  4k  (^“ )  > 
62  =  z(C)  +  (1  +  2k)  (t  +  » 


and,  from  Eq.  (24), 


2  o2(0 

V»  =  “V 


=  i  |~(w  -  3v) (1  -  e~V^)  +  2v(l 
v  |_  3(w  -  v) 


-  e-v«(l  -  .■<) 


Thus,  for  a  narrow,  collimated,  delta-pulse  source  with  energy  E^, 
we  obtain  from  Eq.  (51)  the  irradiance  on  the  exit  plane  z  =  T  as 


J(r,  T;  C)  =  E0  e"^  pr(r,  C)  A(t,  ?)  ,  (52) 

2  2  h 

where  the  radial  distance  r  =  (x  +  y  )  .  Since  x  and  y  are  Gaussian 
random  variables,  r  is  Rayleigh-distributed,  so 


P  (r,  C)  =  —  exp  — § - 

"  V5)  L2ax(° 


COMPARISON  WITH  MONTE  CARLO  SIMULATIONS 

Now  we  digress  to  test  the  validity  of  Eq.  (52).  Consider  the 
power  into  an  infinite-plane,  ir/2-field-of-view  receiver.  For  E^  =  1 
and  K  =  0,  the  received  power  pulse  is  simply  A(T,  £)  since,  from 
Eq.  (52), 


CO  00 

/  J<r'  T;  °  dr  *  A<T>  0  /  Pr(r, 


C)  dr  =  A(t,  £)  . 


*i'*  4  *  afe-  *w  *  ' 
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Bucher  [4]  considers  this  problem  in  his  Monte  Carlo  simulations. 
Figures  10  and  11  plot  the  simulated  and  calculated  pulses  for  T  =  30 
and  80,  respectively.  Although  the  pulse  amplitudes  have  been  normal¬ 
ized  to  unity  in  these  figures,  the  integral  of  the  pulses,  or  received 
energy, 


00 

=  f  a(t, 

r\ 


K)  dC  , 


also  agrees  well  with  Bucher's  Monte  Carlo  curve  fit,  as  shown  in 
Fig.  12. 

Bucher  also  plots  (r)  ,  the  mean  radial  distance  of  exiting 
photons  for  a  cloud  of  optical  thickness  T.  We  calculate 


OO 

/ 


r  hT(r)  dr  , 


(54) 


where  h^r)  dr  is  the  conditional  probability  of  a  photon  exiting 
through  an  annulus  dr  at  r  given  that  it  attains  z  =  T.  Defining  B 
as  the  event  corresponding  to  the  photon  leaving  at  all,  and  using 
conditional  probabilities,  we  have 


h^.(r)  dr  =  ^{r€dr/B}  =  ^{rGdr,  B 


00 
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5°{r€dr,  ?€d?,  B}/f{B} 
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Fig.  11— Comparison  of  simulated  and  calcul 
power  pulses  for  t  =  80 
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But,  from  Eq.  (53),  the  radial  density  of  photons  at  time  ?  is 

2 


^{r€dr,  B/£>)  =  — ^ —  exp 


and 


?{B}  =  f  A(t ,  ?)  d?  =  E, 

J  r\ 


After  substitution,  Eq.  (55)  becomes 


hT(r)  dr  =  Et 
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1  f  x 
Jn  of (5) 


exp 


0  x 
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a(t,  ?)  d? 


and  Eq.  (54)  is  then 
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Reversing  the  order  of  integration  and  recognizing  the  Maxwellian 
density  results  in 


<r)T  =  E"1  /v/2 

=  E"1  /¥?2 


<30 


/  00 
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I  /  r 
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dr 

<  I  o  exp 
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I  4c  “s® 

L  2v^J 
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(?)  A(T,  O  d? 
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t56) 


which  is  numerically  integrated  and  plotted  with  the  Bucher  curve  fit 
in  Fig.  13. 


Fig.  13— Comparison  of  simulated  and  calculated  mean  radial  distance  of 
emerging  photons  versus  optical  thickness  for  (cos  0)  =  0.827 
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A  key  factor  in  optical  communications,  pulse  stretching  is  dra¬ 
matically  affected  by  multipath  delays  due  to  multiple  scattering  in 
clouds,  as  shown  in  Figs.  10  and  11.  To  characterize  multipath  delay 
we  again  adopt  the  mean  delay  =  £  -  T  (discussed  in  Sec.  Ill,  p.  29). 
In  the  diffusion  approximation  for  an  infinite  plane  receiver  at  z  =  T, 

f  ?A(r,  Z)  dZ 

M  -  —Sr - T  •  (57) 

I  A(t,  Z)  dZ 
JZ=o 


Figure  14  plots  against  t  using  Eq.  (57)  and  the  Bucher  simulation 
curve.  The  calculated  pulse  broadens  as  the  diffusion  r.pproximation 
breaks  down  below  t  ***  15. 

Together,  Figs.  10  through  14  illustrate  the  remarkable  corre¬ 
spondence  between  the  calculated  and  the  simulated  results.  The  power 
of  an  analytic  solution  of  the  nonstationary  transport  problem  does 
not,  however,  lie  in  duplicating  highly  aggregated  simulation  results. 
Rather,  once  validated,  the  analytic  model  provides  solution;:  to  prob¬ 
lems  that  would  require  excessive  computer  time  to  simulate  because 
of  sparse  histogram  data.  For  example,  if  the  infinite-plane,  ir/2- 
field-of-v.iew  receiver  is  replaced  by  a  finite  receiver  to  si  udy  field- 
of-view  effects,  then  an  analytic  treatment  may  provide  the  only  prac¬ 
tical  way  to  obtain  results.  We  consider  the  problem  of  finite  receivers 
below. 


(60) 
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where 


p  (r,  O  =  — - r  exp 

<£«>  +  °B 


;«> 4 


(61) 


and  A(T,  £)  is  given  by  Eq.  (51).  The  broad-beam  irradiance  equation 
(60)  is  used  below  to  develop  an  expression  for  the  full  radiance  dis¬ 
tribution  at  the  cloud  exit. 

RADIANCE  AND  RECEIVED  POWER  AT  CLOUD  EXIT 

The  angular  distribution  of  particles  emitted  from  the  exit  sur¬ 
face  of  a  plane-parallel  scattering  medium  with  a  stationary  plane-wave 
input  (Milne's  problem)  is  closely  approximated  by  the  normally 'i 
expression 


2  +  3  cos  0 

g<“> - n - 


0  S  0  S  tt/2,  0  £  (t>  £  2tt  , 


which  is  derived  by  Morse  and  Feshbach  [11]  using  a  steady-state  dif¬ 
fusion  approximation.  Simulations  for  delta-pulse  inputs  have  shown 
that  the  light  emerging  from  an  optically  thick  medium  is  consistent 
with  Eq.  (62)  and  nearly  independent  of  r,  x,  and  £  [4,15].  Thus,  we 
adopt  Eq.  (62)  for  the  angular  distribution  of  photons  leaving  the 
bottom  of  the  cloud,  with  the  source  function  given  by  Eq.  (58). 

The  separability  of  the  spatial  and  angular  dependence  for  light 
emerging  from  the  cloud  bottom  implies  a  radiance  of  the  form 
I(r,  T,  a;  £)  ~  J(r,  x;  £)  g(a).  But,  by  definition, 


J(r,  x;  £) 


=  J  I(r,  x,  as 


£)  cos  9  dco 


where  du)^  =  sin  0  d0  dij)  is  the  solid-angle  element  in  the  direction 
a  =  (0,  <}>) ,  and  the  integral  is  over  ft  =  4 v  steradians.  With  the 
initial  condition  given  by  Eq.  (58),  therefore,  the  radiance  at  the 
cloud  exit  from  Eqs.  (60),  (62),  and  (63)  is 


vklk,  *5*'* 
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I(r,  T,  a;  O  = 


g(a)  cos  9  dw 


^  e"K?(2  +  3  cos  9)  pr(r,  C) 


X  A(t,  o  ,  0  <;  0  £  ir/2  , 


where  Pr(r,  E,)  and  A(t,  £)  are  given  by  Eqs.  (61)  and  (51),  respectively. 

The  power  into  an  on-axis  receiver  at  the  cloud  exit,  with  normal¬ 
ized  aperture  radius  a  and  f ield-of-view  half-angle  9fov>  is 


fov  2tt 


PT(S)  - 


/  /  /  1 

*r=0  0=0  *'d)=0 


cos  9  sin  0  d0  d<J> 


—r-  e  K^(3  -  cos  20,  -  2  cos^  0,  ) 

4  N  fov  fov 


x  I  1  —  exp  {  - 


2  f<4«)  +  oj 


A(T,  O  . 


For  a  =  °°,  0fov  =  tt/2,  Eq  =  1,  and  ic  =  0,  Eq.  (65)  reduces  to 
Pt(£)  =  A(t,  E,) ,  consistent  with  the  infinite-plane-receiver  analysis 
in  Sec.  IV,  p.  43. 


PROPAGATION  BELOW  CLOUD 

Consider  the  narrow  delta  pulse  of  light  emerging  from  the  cloud 

T 

bottom  (z  =  T)  at  position  r*  =  (x',  y')  ,  time  £',  and  direction 
a'  =  (0',  <j>’).  At  time  £  =  K'  +  d  sec  0,  the  pulse  will  illuminate 
the  plane  z  =  T  +  d  at  the  position  _r  =  £*  +  d  tan  0  n  while  main- 
taining  the  direction  a  ~  a',  where  n^  =  (cos  4*’ ,  sin  4* * )  •  Hence 
the  atmospheric  propagator  from  the  plane  z  =  T  to  the  plane  z  =  T  +  d 


is 
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For  calculating  the  power  into  an  on-axis  receiver  at  z  =  T  +  d 
with  aperture  a  and  field-of-view  0£ov»  we  assume  that  the  radiance 
over  the  aperture  is  approximately  equal  to  its  value  at  jr  =  0  when 


a  «  [o2x(Z)  +  a\]h 


for  all  £  >  t.  Thus,  from  Eq.  (68),  the  received  power  is  given  by 
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tan  0 _ 

d  sec  0)  +  a2j 


(2+3  cos  0)  cos  0 


sin  0 
(69) 


which  may  be  evaluated  numerically.  However,  if  0fov  is  also  small, 
Eq.  (69)  reduces  to 


P0  -k£ 

- f  e  A(T,  £  -  d) 

4d^ 
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fov _ 

-  d> +  -2]. 


As  a  check  for  Eq.  (70),  we  obtain 


(70) 
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which  is  identical  to  given  by  Eq.  (65)  evaluated  with  the  usual 
approximations  when  a  and  0£ov  are  assumed  to  be  small. 
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Cry*  *’s ''r.^»S^>‘~- 


To  show  the  effect  of  replacing  Eq.  (69),  which  requires  a  numer¬ 
ical  integration,  with  the  closed-form  equation  (70),  Fig.  15  plots 


ET+d(a-  0fOV>  ' 


PT+d«>  « 


against  0  using  pt+cj(£)  from  both  Eqs.  (69)  and  (70).  The  approxi¬ 
mate  equation  (70)  results  in  under  a  5  percent  error  in  E^^a,  0^^) 
for  0£ov  ^  15  deg  when  Eq  =  1  J,  K=0,  (cos  0)  =  (cos^  0)  =  0.850, 

T  =  30,  d  =  150,  a  =  500,  and  a  =  0.01.  Therefore,  we  recommend  the 
B 

use  of  the  computationally  efficient  Eq.  (70)  whenever  0£ov  ^  15  deg. 
Furthermore,  Fig.  15  indicates  that  beyond  about  8fov  =  20  deg,  the 
signal  energy  into  the  receiver  does  not  increase  appreciably. 

Figures  16  through  21  show  the  effect  of  the  cloud  and  atmospheric 
layers  on  energy  and  pulse  stretching  for  a  finite  receiver  located 
on-axis  a  distance  d  below  the  cloud.  In  each  case,  the  received  energy 
is  given  by 


E  j, (a,  0,  ) 

T+d  f ov 


OO 

■/ 


rT+dro  d5  • 


where  ^T+Cj(?)  is  computed  from  Eq.  (69)  if  0£ov  >  15  deg  or  from  Eq. 
(70)  if  0 £ 0 v  ^  15  deg.  Pulse  stretching  is  again  represented  by  the 
average  multipath  delay 


EWd<a'  9fOV> 


5  PT=d«>  « 


-  (T  +  d) 


We  do  not  parameterize  any  calculations  on  receiver  aperture 
because  of  the  simple  dependence  of  P  (?)  on  aperture  a,  as  evinced 
in  Eqs.  (69)  and  (70).  The  following  input  values  are  used  in  the 

computations  for  Figs.  16  through  21:  E  =  1  J,  ic  =  0,  (cos  0)  = 

2  ^ 

(cos  0)  =  0.850,  and  a_  =  25.  The  thickness  of  the  atmospheric  layer 
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is  expressed  in  kilometers  under  the  assumption  of  a  mean  free  scatter¬ 
ing  length  kg^=  20  m;  in  this  case,  the  standard  deviation  of  the 
incident  beam  is  500  m. 

Figures  16  and  17  show  the  interaction  of  atmospheric  layer  thick¬ 
ness  and  field  of  view.  For  thinner  layers,  wider  fields  of  view  in¬ 
crease  both  received  energy  and  pulse  stretching,  but  the  latter  effect 
is  relatively  slight.  For  thicker  atmospheric  layers  (D  >  3  km),  in¬ 
creasing  field  of  view  beyond  20  or  30  deg  does  not  influence  received 
energy  or  pulse  stretching — because,  if  the  spot  subtended  by  the  re¬ 
ceiver  on  the  cloud  bottom  is  larger  than  the  main  body  of  the  exiting 
beam,  then  increasing  field  of  view  has  little  effect  on  the  received 
pulse.  Figures  18  through  21  display  the  relative  importance  of  cloud 
as  opposed  to  atmospheric  layer  thickness  when  the  receiver  field  of 
view  is  fixed  at  15  deg.  With  regard  to  received  energy  and  especially 
multipath  pulse  stretching,  the  most  important  factor  is  the  optical 
thickness  of  the  cloud,  not  its  distance  above  the  receiver. 
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19— Multipath  time  delay 
15°  field  of  view  and 


versus  distance  below  cloud  for 
various  cloud  thicknesses 
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VI.  IMPLICATIONS  FOR  BLUE-GREEN  PROGRAM 


It  is  difficult  to  assess  the  feasibility  of  satellite-to-mbmarine 
optical  communications  without  a  well-founded,  tractable,  analytic  treat¬ 
ment  for  multiple  scattering  in  clouds.  Simulation  results  have  been 
computationally  expensive,  highly  aggregated  (e.g.,  based  on  infinite 
receiving  planes),  and  often  presented  without  error  bounds.  Though 
for  finite  receivers,  analyses  involving  the  small-angle  assumption  are 
fundamentally  problematical,  especially  when  applied  to  all  but  the 
thinnest  clouds. 

The  mathematical  treatment  of  multiple  scattering  presented  here 
begins  by  deriving  the  scattering  moments — predicated  on  basic  physical 
concepts  and  performed  without  approximations — and  continues  by  develop¬ 
ing  the  diffusion  model,  which  was  validated  by  published  simulation 
results  for  the  identical  problems.  The  credibility  of  both  the  ana¬ 
lytic  model  and  the  simulation  is  established  by  their  remarkable  agree¬ 
ment  for  optical  thicknesses  greater  than  15 — the  range  for  which  clouds 
most  severely  degrade  the  propagation  path.  However,  the  diffusion 
model  provides  expressions  that  are  both  computationally  efficient  and 
cover  a  broad  range  of  problems  at  all  levels  of  aggregation.  The 
moment  formulae  and  the  diffusion  expressions  also  provide  the  physical 
insight  into  the  scattering  process  that  may  be  required  to  resolve  new 
problems  and  existing  anomalies  in  certain  experimental  results,  as 
discussed  below. 

Since  the  diffusion  model  provides  an  expression  for  the  full 
radiance  distribution  on  a  plane  just  above  the  ocean  surface,  it  can 
be  linked  with  other  models  for  radiance  propagation  through  the  ocean 
surface  down  to  a  submarine-mounted  receiver.  With  the  addition  of 
noise,  the  full  propagation  path  would  be  modeled,  allowing  efficient 
t  optimization  of  receiver  aperture  and  field  of  view.  In  any  case,  we 

are  now  prepared  to  comment  on  relevant  parts  of  the  current  model  of 
optical  pulse  propagation  from  satellite  to  submarine. 


NAVY  MODEL 


The  Navy's  blue-green  single-pulse  downlink  propagation  model  [9] 
provides  a  common  base  for  analyzing  satellite -to-submarine  communica¬ 
tion  system  performance.  The  Navy  model  is  modularized  to  permit 
convenient  incorporation  of  new  results  for  each  part  of  the  propaga¬ 
tion  path  (e.g. ,  clouds,  atmosphere,  air/water  interface).  The  de¬ 
velopments  presented  here  are  relevant  to  pulse  shape  and  width,  total 
transmission  through  clouds  and  the  cloud-to-ocean  atmospheric  layer, 
and  the  angular  distribution  of  exiting  radiation  below  the  cloud. 

In  the  following  discussion,  we  either  confirm  the  expressions  used  in 
the  Navy  model  or  suggest  alternatives. 

The  Navy  model  represents  energy  transmission  through  the  cloud 
layer  by  the  Bucher  simulation  fit  [4],  which  is  given  by 


E 


T 


1.69 

vr  +  1.42 


(71) 


in  our  notation.  Figure  12  plots  Eq.  (71)  together  with  the  diffusion 
model  result.  The  close  agreement  between  the  two  curves  demonstrates 
that  Eq.  (71)  is  essentially  correct  and  not  subject  to  the  factor-of- 
two  error  ir.  the  1.69  and  1.42  constants  suggested  in  Ref.  9.  The 
extrapolation  of  Eq.  (71)  to  E  =  1  for  T  =  0,  as  done  in  the  Navy 
model,  seems  appropriate;  moreover,  we  agree  with  the  zenith-angle 
dependence  on  signal  transmission  taken  from  Bucher. 

The  Navy  model  mechanism  for  cloud-to-water  energy  transmission 
appears  to  be  inconsistent  with  the  definition  of  energy  transmission 
through  the  cloud  embodied  in  Eq.  (71).  Indeed,  energy  transmission 
through  the  cloud  can  be  defined  as  the  probability  that  a  photon 
entering  the  cloud  from  the  top  will  eventually  emerge  from  the  bottom. 
Applied  to  energy  transmission  through  the  cloud-to-water  atmospheric 
layer,  this  definition  leads  to  the  conclusion  that  transmission  loss 
is  due  only  to  the  angular  distribution  of  exiting  photons  and  the 
atmospheric  absorption  coefficient.  Therefore,  the  diffusion  model 
calculates  the  transmission  loss  as 
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g(a)  d03a  =  j-  [4  +  ted  (4  -  3i<d)  Ei(-Kd) 


+  3e"Kd(l  -  Kd)]  , 


(72) 


where  the  angular  distribution  g(a)  is  given  by  Eq.  (62)  and  Ei(x)  is 

the  exponential  integral  function  [16].  Equation  (S-6a)  in  Ref.  9 

does  not  account  for  atmospheric  absorption,  is  apparently  based  on 

a  spot  spreading  phenomenon,  assumes  a  Lambertian  distribution  of 

emerging  light,  and  is  substantially  different  from  Eq.  (72). 

■“let 

The  received  pulse  shape  given  by  F(t)  =  t  e  in  Ref.  9  con¬ 
forms  well  with  those  computed  by  Bucher  [4]  and  hence  with  those  com¬ 
puted  here.  For  multipath  time  spreading  due  to  clouds,  we  suggest 
use  of  the  Bucher  simulation  fit 
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1.94 


9 


which  agrees  with  our  results  for  all  T  (Figs.  7  and  14).  Equation  (34), 
which  is  used  in  Ref.  9,  is  not  recommended  because  it  may  not  conform 
to  physical  prerequisites,  as  discussed  in  Sec.  III.  Figures  19  and 
21  show  that  pulse  stretching  due  to  the  cloud-to-water  atmospheric 
path  is  negligible  compared  with  that  due  to  clouds.  Hence  we  agree 
with  Eq.  (S-25b),  which  sets  cloud-to-water  pulse  stretching  to  zero. 


EXPERIMENTAL  VERIFICATION 

Because  the  diffusion  model  provides  expressions  for  the  power 
received  by  a  finite  aperture  and  f ield-of-view  receiver  [Eqs.  (69) 
and  (70)]  located  an  arbitrary  distance  below  a  given  cloud  layer,  the 
model  is  suitable  for  verification  by  comparison  with  experimental 
results.  Indeed,  an  experiment  is  currently  under  way  on  the  island 
of  Kauai  partly  for  this  purpose.  However,  experimental  verification 
is  fraught  with  complications,  such  as  (1)  the  difficulty  of  measuring 
experimental  input  parameters  (e.g.,  T,  (cos  6) ,  k  )  accurately  enough 

cL 

to  conclude  that  the  model  is  correct,  and  (2)  the  inconsistency  of 
the  geometric  boundaries  governing  the  experiment  with  the  assumptions 
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used  to  develop  the  model  (e.g.,  a  plane-parallel  cloud  layer).  Ob¬ 
taining  accurate  measurements  of  physical  parameters  is  beyond  the 
scope  of  this  report,  but  the  effect  of  boundary  problems  is  discussed 
below. 

One  important,  recent  multiple  scattering  experiment  was  the  Naval 
Ocean  Systems  Center  (NOSC)  work  on  fog  propagation  reported  in  Ref.  1. 
The  effect  of  the  ground  as  an  absorbing  boundary  for  photons  has  been 
of  some  concern  in  interpreting  the  results  of  the  NOSC  experiment. 

By  using  the  boundary  value  techniques  discussed  in  Sec.  IV,  the  dif¬ 
fusion  model  can  be  extended  to  account  for  an  absorbing  plane  parallel 
to  and  laterally  displaced  from  the  axis  of  propagation.  Similar  ex¬ 
tensions  can  be  developed  to  model  rectangular  clouds,  or  rectangular 
slots  or  holes  in  clouds. 

A  scaled-down  multiple  scattering  laboratory  experiment  has  been 
suggested,  using  large  tanks  filled  with  water  and  suspended  oil 
droplets.  The  dimensions  of  the  tank  and  the  absorptive  properties 
of  the  walls  required  to  model  a  cloud  without  horizontal  boundaries 
are  two  of  the  technical  issues  that  should  be  resolved  before  such 
an  experiment  is  undertaken.  Our  results  indicate  that  tank  boundaries 
(perpendicular  to  the  coordinate  frame  of  photon  motion)  can  be  modeled 
in  the  diffusion  theory  sense,  as  either  absorbing,  reflecting,  or 
elastic,  so  long  as  photons  enter  the  tank  at  least  15  optical  thick¬ 
nesses  from  the  side  walls.  Thus,  the  tank  dimensions  and  side-wall 
absorptive  properties  may  not  inhibit  the  laboiatory  verification  of 
the  diffusion  model. 


J 


1 1  -  f»  •* *jf«WAV '  -  Wp ~  aOP'V’J  <*J  *  W  • 


-71- 


REFERENCES 


»  v> 


£ 


$ 


I 


1.  Mooradian,  G.  C.,  et  al.,  "Blue-Green  Pulsed  Propagation  through 

Fog,"  Applied  Optics,  Vol.  18,  No.  4,  15  February  1979,  pp. 
429-441. 

2.  Chandrasekhar,  S.,  Radiative  Transfer,  Oxford  at  the  Clarendon 

Press,  Oxford,  i  s  'nd,  1950. 

3.  Stotts,  L.  B.,  "The  Radiance  Produced  by  Laser  Radiation  Trans- 

versing  a  Particulate  Multiple-Scattering  Medium,"  J.  Opt.  Soc. 
Am.,  Vol.  67,  No.  6,  June  1977,  pp.  815-819, 

4.  Bucher,  E.  A.,  "Computer  Simulation  of  Light  Pulse  Propagation 

for  Communication  through  Thick  Clouds,"  Applied  Optics,  Vol. 

12,  No.  10,  October  1973,  pp.  2391-2400. 

5.  Parzen,  E. ,  Stochastic  Processes,  Holden-Day,  San  Francisco,  1962. 

6.  Reif,  F.,  Statistical  Thermal  Physics,  McGraw-Hill,  New  York,  1965. 

7.  Arnush,  D.,  "Underwater  Light-Beam  Propagation  in  the  Small-Angle- 

Scattering  Approximation."  J.  Opt.  Soc.  Am.,  Vol.  62,  No.  9, 
September  1972,  pp.  1109-1111. 

8.  Stotts,  L.  B.,  "Closed  Form  Expression  for  Optical  Pulse  Broaden¬ 

ing  in  Multiple-Scattering  Media,"  Applied  Optics,  Vol.  17, 

No.  4,  15  February  1978,  pp.  504-505. 

9.  Technical  Advisor  to  the  Blue-Green  Optical  Communication  Pro¬ 

gram  Joint  Coordinating  Committee,  Naval  Blue-Green  Single- 
Pulse  Downlink  Propagation  Model,  Naval  Ocean  Systems  Center, 

San  Diego,  Calif.,  TR  387,  1  January  1979. 

10.  Feller,  W. ,  An  Introduction  to  Probability  Theory  and  Its  Appli¬ 

cation/,  Vol.  2,  Wiley,  New  York,  1966. 

11.  Morse,  P. ,  and  H.  Feshbach,  Methods  of  Theoretical  Physics,  McGraw- 

Hill,  New  York,  1953. 

12.  DeGroot,  J.  S.,  R.  M.  Frank,  and  H.  ,S.  Stewart,  Effects  of  Clouds 

on  X-Ray  Induced  Fluorescent  Light  from  a  Nuclear  Burst  in  Space, 
E.  H.  Plesset  Associates,  Santa  Monica,  Calif.,  C73-65(u)ll, 

29  January  1965. 

13.  Fritz,  S.,  "Scattering  and  Absorption  of  Solar  Energy  by  Clouds," 

Ph.D.  Thesis,  Massachusetts  Institute  of  Technology,  Cambridge, 
1953. 


I 


PrfLCLuL.G  br.C£.  bL«i<K-N0T  FlilffiD 


j! 


14.  Carslaw,  H.  S.,  and  J.  C.  Jaeger,  Conduction  of  Heat  in  Solids 3 

2d  ed.,  Oxford  at  the  Clarendon  Press,  Oxford,  England,  1959. 

15.  Ciany,  C.,  G.  Schroeder,  and  G.  Lee,  "Optimum  System  Design  for 

Optical  Satellite  to  Submarine  Communication  Links,"  McDonnell 
Douglas  Astronautics,  St.  Louis,  Mo.  (unpublished). 

16.  Grads'nteyn,  I.  S.,  and  I.  M.  Ryzhik,  Table  of  Integrals ,  Series s 

and  Products ,  4th  ed.,  prepared  by  Yu.  V.  Geronimus  and  M.  Yu. 
Tseytlin,  translation  edited  by  Alan  Jeffrey,  Academic  Press, 
New  York,  1965. 


